clc
clear

z_max = 40000;
nc    = 40; 
m     = 2;

nl    = nc+1; % number of levels(cell boundaries) on vertical distribution
eta_max = z_max;
deta  = z_max / nc;
eta   = 0+deta/2:deta:eta_max-deta/2;
eta   = eta';

x = ( eta / eta_max - 1 / 2 ) * m;
dxdeta = m / eta_max;

x_min = -0.5*m;
x_max = 0.5*m;
dxh   = x_max - x_min;

coef = 0.5  / x_max * z_max;

f = coef * ( tanh(x) + 1 ) * dxdeta;
df = coef * ( 1 - tanh(x).^2 ) * dxdeta^2;
ddf = coef * ( - 2 * tanh(x) .* ( 1 - tanh(x).^2 ) ) * dxdeta^3;

for i = 1:nc
    F(i,1) = coef * ( ( log(cosh(x(i)) ) - log( cosh(x_min) ) ) + x(i) + dxh / 2 );
end

figure 
hold on
plot(eta,F,'Color','m','Marker','o')
plot(eta,f,'Color','r','Marker','o')
plot(eta,df,'Color','b','Marker','o')
plot(eta,ddf,'Color','k','Marker','o')

level = zeros(2,nl);
for i = 1:nc
    level(:,i) = F(i);
end

figure
plot(level,'b')